A&A manuscript no. 

(will be inserted by hand later) 



Your thesaurus codes are: 
02(02. 01. 1;02.09.1;02. 19.1) 



ASTRONOMY 
AND 
ASTROPHYSICS 
1.2.2008 



OS 
Os 



Acoustic emission and corrugational instability of shocks 
modified by strong particle acceleration 



M. MondQ and L. O'C. Drury 

Dublin Institute for Advanced Studies 
5 Merrion Square 
Dublin 2, Ireland 

Received; Accepted 



o 

Q 

00 



> 



(N 

o> 



Of 
i 1 

o 



X 



Abstract. The effect of particles that undergo strong 
diffusive-shock-acceleration on the stability of the acceler- 
ating shock is investigated. A two-fluid model is employed 
in which the accelerated particles are treated as a fluid 
whose effect is incorporated as an additional pressure in 
the momentum equation. The Dyakov and Kontorovich 
stability criteria are used in order to study the stability 
of those steady-state shocks that contain a gas sub-shock. 
The downstream conditions of the latter are parametrized 
by the ratio of the upstream acclerated-particles pressure 
to the total pressure. For some range of values of that 
parameter, three possible downstream states are possible 
for each upstream state. It is shown that in that range of 
parameters the shocks are either corrugationally unstable 
or lose their energy by spontaneous emission of acoustic 
as well as entropy- vortex waves. 



1. Introduction 

Twenty years ago a process, now generally called diffusive 
shock acceleration, was described whereby a significant 
part of the kinetic energy flowing into an astrophysical 
shoc k cou ld be used to accelerate charg ed par ticles (Krym- 
sky, |1977| ; Bell 1978a,b; Axford et al., |l977t Blandford & 
Ostriker, 197S). Clearly the reaction of the accelerated 
particles on the shock dynamics and structure then be- 
comes an important effect. In the intervening years much 
work has been done, but the problem is a difficult one 
and many aspects remain open (see reviews by Jones & 
Ellison, |l99l| ; Bere zhko fc Krymsky, [198$ Blandford & 
Eichler, 1987; Volk, 1987). A significant recent advance is 



the approximate analytic theory of Malkov (1997 
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For understandable reasons most of this work has as- 
sumed steady, or at least quasi-stationary, shock struc- 
tures and there has been relatively little work on the ques- 
tion of shock stability. In this pa per we apply the classi- 
cal cri teria developed by Dyakov ( 1954 ) and Kontorovich 
( 1957 ) for the corrugational instability of shock fronts and 
the spontaneous emission of sound and entropy waves to 
the simplest model of shocks strongly modified by par- 
ticle acceleration, the two-fluid model of Drury & Volk 



(1981). While this model has obvious defects (Heavens & 
Drury, 198$ Achterberg et al., 1984) and has clear limi- 
tations (Drury et al., 1995) it has proven to be a useful 
tool in interpreting numerical studies (Falle & Giddings, 
1987 ; Jones & Kang, 199C ) when used with appropriate 
caution. 



2. The two-fluid model 

The key assumption in the two-fluid model is that the 
spatial transport of the accelerated particles can be rep- 
resented by diffusion with a single effective diffusion co- 
efficient. This enables one to integrate over the particle 
energy spectra and regard the accelerated particles as a 
second fluid with significant energy density and pressure, 
but negligible inertia. For historical reasons the acceler- 
ated particles are often identified as cosmic rays and the 
energy density and pressure denoted Eq and Pq, 



Er I 4np 2 T(p)f(p)dp 



Pi 



c 



Airp 



-f(p)dp 



(1) 
(2) 



where p is the magnitude of the accelerated particle mo- 
mentum, f(p) is the isotropic part of the phase space den- 
sity, T(p) the associated kinetic energy and v(j>) the par- 
ticle speed. 

The two fluid model then consists of the standard hy- 
drodynamical equations for the background gas with den- 
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sity p, velocity U, energy density P G and pressure P G in- 
cluding both pressures in the momentum equation and a 
diffusive energy flux term with coefficient k in the internal 
energy equation for the accelerated particle fluid 
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r\ 1 

-U + U-VU + -V (P c + P G ) 
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P c + U • VP C + 7cPcV • U = V (kVP c ) . 



(3) 
(4) 
(5) 
(6) 



At discontinuities this system has to be replaced by the 
usual conditions of mass, momentum and energy conser- 
vation and continuity of Pc . The adiabatic exponents 7 G 
and 7c are defined by 



Pg,c = ( 7 G,c - 1)Pg,C- 



(7) 



For an ideal monatomic gas 7 G = 5/3 and 4/3 < 7c < 
5/3. 

The simplicity of the two-fluid model is rather decep- 
tive. All information about the energy (or momentum) 
spectrum of the accelerated particles has been put in the 
two closure parameters k and 7c and, perhaps not sur- 
prisingly, the solutions turn out to depend very sensitively 
on the assumptions made for these parameters. To some 
extent this can be overcome by using simple physical mod- 
els for the spectral evolution to improve the estimates 



(Duffy et al., |f994j ). In addition Malkov & Volk ( pL996[ ) 
have shown that even in circumstances where the simple 
two-fluid model is invalid a form of renormalised two-fluid 
model is applicable. 

Our approach in this paper is rather different, We re- 
gard the two-fluid equations not as an approximation to 
some more complicated system, but as the simplest model 
system (in currently popular terminology a "toy" model) 
in which we can study the interaction between particle 
acceleration and shock structure. We note that if the ac- 
celerated particles were to have an energy independent 
diffusion coefficient the two-fluid equations would be ex- 
act. Thus the model is physical and does not violate any 
fundamental physical principles. 

All possible stationary shock structures within this 
model were classified by Drury & Volk, (1981). It is 
interesting to recall one of the concluding remarks in 
this paper. "Perhaps the most interesting, certainly the 
most characteristically nonlinear, feature of this model is 
that under certain circumstances it has three solutions. 
Whether all are stable and can occur as time asymptotic 
states in a physically reasonable evolution remains to be 
studied...". This paper is a belated answer to this ques- 
tion. 



3. Stability criteria for shocks 

The study of the s tability of shocks was first undertaken 
by Dyakov, ( |I954| ). He considered a planar shock wave 
propagating in an unbounded medium and investigated 
the evolution of small sinusoidal corrugations of the shock 
front. Dyakov showed that, under certain conditions, per- 
turbations that are confined to the front's area grow expo- 
nentially with time. This phenomena is termed corruga- 
tional instability. Shocks that are corrugationally unstable 
can only be short-lived because a reorganization of the en- 
tire flow will take place on a time-scale of the order of the 
inverse growth rate. 

In addition to the corrugational instability, Dyakov 



( |1954D and later Kontorovich ( |1957|]1957[ ) determined the 
conditions under which small acoustic as well as entropy- 
vortex perturbations in the form of sinusoidal two di- 
mensional waves can be emitted from the shock front. 
This phenomenon is termed spontaneous emission. Even 
though the spontaneously emitted waves do not grow with 
time, their occurence may eventually render the flow be- 
hind the spontaneously emitting shock unstable; this is 
because the shock's energy is continuously carried away 
by the outgoing waves and this finally results in the reor- 
ganization of the flow. 

Both corrugational instability and spontaneous emis- 
sion may be treated as resonant reflection of acoustic 
waves from the shock front where the reflection coefficient 
of the acoustic waves from the shock is infinite. Wereas 
historically the reflection coefficient has been calculated in 
terms of the geometrical parameters of the problem, such 
as the angles of incidence and reflection, Mond & Rutke- 
vich (1994) have recently recast the traditional Dyakov- 
Kontorovich formulation into a frequency-dependent rep- 
resentation of the reflection coefficient. Such a representa- 
tion is convenient in solving various boundary value prob- 
lems such as, for example, investigation of the eigenper- 
turbations between the shock and some other reflecting 
surface behind it. A further advantage of the frequency- 
dependent formulation is evident when dissipative effects 
are taken into account and the relevant frequencies be- 
come complex. In this case simple geometric concepts like 
the angle of incidence are of restricted physical meaning. 

The reflection coefficient is calculated by considering 
an incident acoustic wave as well as reflected acoustic and 
entropy-vortex waves in the region downstream from the 
shock of the form 



Sf exp(—iwt + ik • r), 



(8) 



where Sf is the amplitude of the perturbation of the rel- 
evant physical quantity. While u> and the component of 
k which is parallel to the unperturbed shock front are 
the same for all the waves involved in the process, the 
perpendicular component of k is determined seperately 
for each wave according to the appropriate dispersion re- 
lation. The upstream perturbations are zero due to the 



M. Mond & L. O'C. Drury 



3 



supersonic velocity there. The reflection coefficient is ob- 
tained by imposing the continuity of the tangential veloc- 
ity component as well as the appropriate jump condition 
of the normal velocity components at the perturbed shock 
front. Its frequency-dependent representation is given by: 



where 



K 



f(g) (g 2 - l) 1 / 2 
f(q) + («Z 2 - 1)V2 



where 



(1 + fo)r?M 2 
2g(l-M 2 ) : 



q = 



k\\c 2 (l 



Mi) 1 / 2 



(9) 



,(10) 



Mi is the downstream Mach number and u> and fcii are 
the frequency, and the wave vector component parallel to 
the shock's front, of the acoustic as well as of the entropy- 
vortex waves. As will be seen later, the parameter h plays 
a crucial role in determining the shock's stability and is 
given by 



h = 



dp 2 
dp2 



(11) 



where v, p and p are the velocity, density and pressure of 
the gas, respectively and the subscripts 1 (2) denotes the 
corresponding values in the upstream (downstream) side 
of the shock. Thus, h is given in terms of the derivative of 
the downstream density with respect to the downstream 
pressure for fixed upstream conditions, i.e., in terms of the 
slope of the Hugoniot curve. 

The problem of stability may be cast now into the fol- 
lowing question: under what conditions does the denomi- 
nator on the right hand side of eq. (^) become zero? The 
answer is that two families of two dimensional acoustic 
as well as entropy-vortex waves may give rise to an in- 
finite reflection coefficient under two distinct conditions. 
The first family of waves is characterized by purely imag- 
inary ui and purely imaginary k± (the wave vector com- 
ponent perpendicular to the shock's front). Hence, such 
waves grow exponentially in time and decay exponentially 
away from the shock. They give rise to the corrugational 
instability and the conditions for their existence are 



h < -1, 



h>l + 2M 2 



(12) 



When condition (113) is satisfied, the requirement of an 
infinite reflection coefficient results in a linear dispersion 
relation between the growth rate of the instability and the 
corrugation wavelength, kn . 

The second family of waves that can give rise to an 
infinite reflection coefficient is characterized by real fre- 
quencies and real wave vectors. Those waves are the man- 
ifestation of spontaneous emission and the condition for 
their occurence is 



kr = 



l-(l + r?)M 2 
1 + (77 - l)Af| 



(14) 



where 77 = p 2 / Pi is the compression ratio. The implication 
of condition ( |l~3| ) is that if it is satisfied, there is a one 
parameter family containing an infinite number of two di- 
mensional acoustic as well as entropy- vortex waves any of 
which can be spontaneously emitted from the shock front. 
A linear dispersion relation exists between the frequen- 
cies of the spontaneously emitted waves and their parallel 
wave vectors. 

4. Stability of shocks modified by particle acceler- 
ation 



h c < h < 1 + 2M 2 



(13) 



When conditions (|l2j) and (13) are applied to shocks in an 
ideal polytropic gas the result is absolute stability against 
both corrugational instability and spontaneous wave emis- 
sion. This result is not unexpected and has been well 
known experimentally since the early days of shock wave 
research. However, as nonideal processes become progres- 
sively more important, the stability properties of the flow 
may significantly change and the shock may become sus- 
ceptible to one of the instabilities discussed in the pre- 
vious section. An example for such occu rence has been 
recently discussed in Mond & Rutkevich ( |l997|) where it 
has been shown that strong ionizing shocks spontaneously 
emit acoustic waves if their Mach number exceeds a cer- 
tain value. 

Here, the effect of shock accelerated particles on the 
stability of the accelerating shock is investigated. For that 
purpose, the two-fluid model discussed in section 2 is em- 
ployed. Within the framework of this model, a shock is a 
transition layer between two uniform states, whose length 
is determined by the mean diffusion coefficient of the par- 
ticles. Following Drury (1983), it can be shown that the 
transition layer may either be smooth or contain a sub- 
shock in the background gas. It is the latter case that will 
be the focus of the current investigation. 

For the purpose of this stability analysis the transition 
layer that includes the gas sub-shock will be regarded as a 
single surface of discontinuity. It is this surface that will be 
called the shock from now on. Such a procedure is admis- 
sible if the wavelengths of the investigated perturbations 
are much larger than the thickness of the transition layer. 

4-1- The Hugoniot curve 

It is obvious from equations (pT|)-(p^|) that the shape of 
the Hugoniot curve plays a crucial role in determining 
the stability properties of shocks that propagate into a 
given uniform state. In order to obtain the Hugoniot curve, 
the conservation equations that relate the flow variables 
on both uniform-state-sides of the shock are written in a 
frame of reference that is moving with the shock's velocity. 
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For that purpose it is convenient to introduce the following 
nondimensional variables: 



El 
Pi 



n=-—. £ 

and 

A G = 



Pci 



Pgi + Pci : 



P G 2 



Pgi + Pci ' 



Ac = 



P C 2 



Pgi + Pci 



(15) 



(16) 



where the subscript 1 (2) denotes the values of the up- 
stream (downstream) corresponding variables. Setting all 
time derivatives to zero in equations (3)-(6), the following 
equations are obtained: 





J 2 (l--)-A G -A c = 
V 


-i, 


(17) 




1 x r A G A c 

- —z) - 1 G 1 C — 


-o, 


(18) 




J 2 (---) + A a -A G = 

Va V 


0, 


(19) 


-J 2 (- 




0, 


(20) 




A a = 


v2 G - 


(21) 


where 9 


= (r G (i-o+r c O,r GlC 


= lG,c/{lG,C- 


1) and 



J is the normalized mass flux and is given by 

Pgi + Pci' K ' 

Equations (17) and (18) represent the momentum and 
energy conservation, respectively, across the shock while 
equations (19) and (20) represent the conservation of the 
background gas momentum and energy across the gas sub- 
shock. The variables A a and rj a are the normalized pres- 
sure and density at the foot of the gas sub-shock and are 
adiabaticaly related to the upstream conditions according 
to eq. (21). 

The solution of the system of equations (17)-(21) re- 
sult in a one parameter (£) family of Hugoniot curves. Two 
such typical curves for 7 G = 5/3 and 7c =4/3 are shown 
in Fig. 1 for £ = 0.45 and in Fig. 2 for £ = 0.1. As can be 
seen in both figures, for small values of Mi the Hugoniot 
curve follows its single-fluid, 7 = 5/3 counterpart while 
for large values of Mi it is cosmic-rays dominated as it 
asymptotically approaches the value of 7 as predicted for 
a single-fluid with 7 = 4/3. It should be noted, however, 
that for high enough values of Mi no physically accept- 
able solution of the system of equations (17)- (21) exists 
which indicates that only a smooth transition layer (with- 
out a gas sub-shock) may exist for that range of parame- 
ters. The two asymptotic parts of the Hugoniot curve are 
connected by an intermediate section. As can be seen in 
Fig. 2, for small values of £ for some portion of the in- 
termediated section there are three possible downstream 
states for each upstream state. As will be seen later on, 
the shock instabilities occur at that intermediate range of 
parameters. 



4 5 6 7 8 9 

Fig. 1. The Hugoniot curve for £ = 0.45 
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Fig. 2. The Hugoniot curve for £ = 0.1 



4-2. The stability analysis 

As was discussed in the previous section the concept of 
a single surface of discontinuity is meaningful only for 
waves whose wavelengths are much longer than the thick- 
ness of the transition layer. In that long wavelength limit 
the acoustic perturbations in the background gas couple 
to the accelerated particles and travel at the enhanced 
speed 



c 2 = 



JgPg + IcPc 



(23) 



In this case the parameters h and M2 that are needed 



for the calculation of the stability criteria (12) and (13) 
may be expressed in terms of the nondimensional variables 
defined in Eqs. (fjJI) and dlrf ) as 



h = - 



J 2 dr/ 
Y<9A 



(24) 
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where A = Aq + Ac is the normalized total downstream 
pressure, and 



J 



77(r G A G + r c A c ) 



(25) 



The variables h and M2 were calculated according to 
( pij ) and(^5|) along the Hugoniot curves and then were 
used in criteria (|l^) and ([l3|). The results are shown in 
Fig. 3 and 4 where h c — h is plotted for £ = 0.45 and 
£ = 0.1, respectively. 




Fig. 3. The behaviour of h c — h as a function of 77 for £ = 0.45 
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Fig. 4. The behaviour of h c — h as a function of r\ for f = 0.1 



It is obvious that all possible shocks for the £ = 0.45 
case are stable against both corrugation instability as 
well as spontaneous emission. This changes however for 
£ = 0.1. In that case there are three possible downstream 



states for each upstream state at some intermediate sec- 
tion of the Hugoniot curve. According to Fig. 4 shocks that 
belong to the descending part of the intermediate section 
are unstable under spontaneous acoustic emission. This 
is so since along that portion of the Hugoniot curve h is 
positive (see Eq. (Q)) while h c remains negative. Further- 
more, Fig. 2 indicates that the shocks near the "knees" of 
the Hugoniot curve are corrugationally unstable. An al- 
ternative algebraic form for the parameter h is easily seen 
to be 



h = -l/[l + 2(r] - 1) 



7? dUi 



(26) 



This shows that that h = — 1, indicating the onset of cor- 
rugational instability, at precisely those points where the 
graph of 77 against U\ is vertical. A slight further inflec- 
tion of the graph then allows the denominator to become 
zero and h — > 00. For such high values of h the dispersion 
relation that is obtained for the corrugational instability 
is approximately given by 



1 



(27) 



so that the growth rate according to the definition in eq. 

Tol) is 



(28) 



a = ^r}M 2 k\\C2. 



To complete the stability analysis it should be noted 
that in the short wavelength limit the acoustic waves de- 
couple from the particles and propagate within the back- 
ground gas at the gas sound speed (jcPg/p) 1 ^ 2 ■ In this 
limit the acoustic waves propagate through the transition 
layer and are reflected off the gas sub-shock. Hence, the 
flow is stable under short wavelength perturbations. 

5. The effect of particle diffusion 

As was discussed in section 2 the spontaneously emitted 
waves are marginally stable in the sense that their fre- 
quencies are purely real. It is of interest, therefore, to 
investigate the effects of dissipative mechanisms on the 
stability of the flow as they give rise to imaginary parts 
of both the frequencies and the perpendicular wave vec- 
tor. Hence the question to be asked is whether, under 
the effects of dissipation, some of the spontaneously emit- 
ted eigenperturbations acquire positive imaginary parts of 
both their eigenfrequencies and their perpendicular wave 
vectors. This question is still an open one and no general 
satisfactory answers are available. Furthermore, conflict- 
ing claims may be found in the literature with respect to 
the stabilizing/destabilizing effects of dissipation. 

In order to study the effect of the particles dissipation 
let us reconsider the denominator of the right hand side of 
Eq. (J^j. It was tacitly assumed before that the coefficients 
in that expression (and in particular Mi) were real. Then 
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the resulting cigcnfrequencies were either purely imagi- 
nary (corrugation instability) or purely real (spontaneous 
emission). However, when the effects of particle diffusion 
are included this ceases to be true and the Mach num- 
ber, and consequently the frequencies of the spontaneously 
emitted waves, become complex. Linearizing Eqs. (3)-(7), 
assuming that the diffusion is small and that kn >> k± 
results in the following expression for the sound speed 



IcPc 



IgPg + IcPc 



(29) 



where cq is the sound speed in the absence of diffusion and 
is given by (E3h and x is defined by 



Km K 
COI 



(30) 



With the above modified expression of the sound speed 
the Mach number is given by 



M 2 = M 2 (l + i 



IcPc 



IgPg + IcPc 



x). 



(31) 



Using expression ( |3l| ) the denominator of the right hand 
side of Eq. (g) is equated now to zero and solutions for 
the frequencies are searched such that the corresponding 
imaginary parts of the perpendicular wave vectors are pos- 
itive. The sign of the imaginary part of the perpendicular 
wave vector is checked by using the dispersion relation for 
the acoustic waves that propagate away from the shock 



w[-M +(l-g- 2 ) 1 / 2 ] 
c(l - M 2 ) 



(32) 



Using the results for a shock (Mi — 4.5) that belongs to 
the spontaneous emission regim it is found that the par- 
ticles dissipation give rise to the damping of the sponta- 
neously emitted waves. As is expected, the damping rate 
is proportional to fcjj and is much smaller that the cor- 
responding real part of the frequency. Thus, the energy 
of the shock is continuously being carried away from the 
shock by the spontaneously emitted waves and is subse- 
quently deposited over a length scale of l/k± . 

6. Discussion 

This paper fills a significant gap in the theoretical under- 
standing of the two fluid model by demonstrating what 
had previously only been conjectured, that in those cases 
where the model produces multiple solutions the interme- 
diate solution is unstable. It is also, as far as we know, 
the only case where Dyakov's criterion for corrugational 
instability is actually realised in a physically plausible sys- 
tem. In addition to its theoretical interest, we note that 
the two fluid model is quite extensively used as the ba- 
sis for numerical studies of particle acceleration effects in 
astrophysical systems; clearly the knowledge that shock 



instabilities of the type described here can, indeed must, 
occur in such simulations is an important piece of infor- 
mation. It will be interesting to see whether similar insta- 
bilties are found in more realistic models with momentum 
dependent diffusion. As we have already noted the two 
fluid model captures much of the correct physics and is 
not unphysical. This, together with Malkov's renormali- 
sation interpretation of the two fluid model, suggests that 
these effects should also occur in the more complicated 
models. 
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